#include "dyn_quad_qSimpson.h"

qSimpson::qSimpson(Functor &_func, double _a, double _b, double _tol, int _maxIterations)
    : Quad(_func, _a,_b, _tol)
{
    setMaxIterations(_maxIterations);
}

double qSimpson::doQuadrature()
{
    // We prepare values to hold the result
    double s;           // Actual result
    double sOld=0.0;    // Previous result.
    double sTemp;       // Temporary result to compute the extrapolation (S_2N)
    double sTempOld=0.0;// Temporary old result.

    // The driver is defined.
    /*! \todo Same as qTrap. */
    TrapzDriver t(func,a,b);

    // We compute the successive refinements until convergence.
    for (int i=0;i<maxIterations;i++)
    {
        // We compute the integral.
        sTemp = t.computeNextRefinement();

        // We extrapolate the results via Simpson's formula.
        s = (4.0*sTemp-sTempOld)/3.0;

        // After a determined number of iterations,
        // we check for convergence (this avoids spurious convergence
        // due to small stepsize).
        if (i>5)
                if (fabs(s-sOld) < tol*fabs(sOld) || (s == 0.0 && sOld == 0.0))
                    return s;

        sOld=s,sTempOld=sTemp;
    }
    cout << "The algorithm has not converged." << endl;
    return 0.0;
}
